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Quantum Treatment for Bose-Einstein Condensation in Non-Equilibrium Systems 
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We develop an approach based on stochastic quantum trajectories for an incoherently pumped 
system of interacting bosons relaxing their energy in a thermal reservoir. Our approach enables the 
study of the versatile coherence properties of the system. We apply the model to exciton polaritons 
in a semiconductor microcavity. Our results demonstrate the onset of macroscopic occupation in 
the lowest-energy mode accompanied by the establishment of both temporal and spatial coherence. 

We show that temporal coherence exhibits a transition from a thermal to coherent statistics and 
the spatial coherence reveals off-diagonal long-range order. 
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Introduction .— A typical signature of Bose-Einstein 
condensation is the formation of macroscopic occupa¬ 
tion in the single-particle ground state of a many-body 
system in thermal equilibrium. This collective state ex¬ 
hibits distinctive spatial and temporal coherence proper¬ 
ties. In solid-state systems, however, the bosons have a 
finite lifetime and hence we need to study steady-state 
properties instead of thermal equilibrium pQ. We re¬ 
fer to this scenario as quasi-Bose-Eistein condensation 
(qBEC). Here, the macroscopic ground-state occupation 
is induced by the relaxation of the higher excited states 
which are pumped by some source reservoir. Such qBEC 
occurs in systems of magnons [2] 13 ] , indirect excitons g], 
or exciton polaritons GHZ] under nonresonant excitation. 

Exciton-polaritons arise from the strong light-matter 
coupling enhanced in semiconductor microcavities [8]. 
They behave like bosonic quasiparticles at moderate con¬ 
centrations (< 10 11 cm -2 ). Due to their small effective 
mass polaritons can manifest quantum coherent prop¬ 
erties up to room temperatures in wide-bandgap mate¬ 
rials ® mg. Polariton qBEC emerges as the result of 
boson-boson interactions and energy exchange with the 
environment. In particular, the scattering of polaritons 
with acoustic phonons m plays a key role, as demon¬ 
strated in a number of recent experiments mm- 

The formation of the quasicondensate is associated 
with emission of coherent laser-like light from the micro¬ 
cavity ini mg. However, such emission is not sufficient 
to prove the existence of qBEC m and since the spatial 
and temporal coherence properties need to be addressed 
in more detail. Experimentally, the temporal coherence is 
described by the second-order temporal coherence func¬ 
tion, g( 2 )(r), where r is the delay between two photode¬ 
tection events. In particular, g^ 2 \r = 0) exhibits a tran¬ 
sition from a thermal g^ 2 \0) = 2, to a coherent statistics, 
g( 2 )(0) = 1 mm- The spatial coherence, also referred 
to as the off-diagonal long-range order [T9U23] , is charac¬ 


terized by a slowly decaying first order spatial coherence 
function, g^\Ax), between regions separated by Ax. 

Whereas experimental techniques based on Michelson 
interferometry and Hanbury Brown and Twiss (HBT) 
setups are well established to measure spatial and tem¬ 
poral coherence, a general theory accounting for many- 
body quantum correlations, environmental interactions, 
and photon counting is lacking. On one hand, the widely 
used semiclassical Boltzmann equation approach [241128] 
to the polariton dynamics is based on the assumption of 
complete incoherence of the system. Thus, the quan¬ 
tum states of the system are taken to be completely 
uncorrelated. On the other hand, approaches based 
on the Gross-Pitaevskii equation [29 under resonant or 
nonresonant excitation [n ed] are successful in explain¬ 
ing a plethora of recent experiments [3TH33] . but they 
assume global coherence and therefore cannot describe 
phonon-assisted relaxation. The truncated Wigner ap- 
proches [34], which involve additional noise terms in the 
Gross-Pitaevskii equation, are based on several limiting 
assumptions and are not designed to describe multimode 
systems. Recently, Boltzmann and Gross-Pitaevskii 
equations have been merged in a classical treatment [35] 
which, unfortunately, does not provide an accurate de¬ 
scription of the onset of coherence. 

Master equation approaches allow to account for both 
coherent and incoherent processes in the polariton dy¬ 
namics [36] and spatial coherence has been recently an¬ 
alyzed in a one-dimensional case gam]. However, such 
model, involving a cumbersome hierarchy of coupled and 
truncated equations, becomes computationally very de¬ 
manding in higher dimensions, and its application seems 
to be restricted to ID structures. 

The model considered in Ref. [39] includes energy re¬ 
laxation in a phenomenological way, operate with a clas¬ 
sical stochastic field, and require unknown fitting param¬ 
eters. Thus they seem unable to describe the desired co- 
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herence properties. Furthermore, the long-range interac¬ 
tions prevent the use of powerful quantum methods based 
on the density matrix renormalization group theory m • 

In this Letter, we develop a highly parallelizable quan¬ 
tum stochastic approach (41 j going far beyond the sin¬ 
gle mode description of Ref. [421 to describe incoherently 
driven interacting bosons with dissipation caused by a 
thermal reservoir. The formalism is based on stochas¬ 
tic evolution of the multimode system wave function in 
its full Hilbert space. It allows the reconstruction of the 
system density matrix from which correlations such as 
<? (2) ( 0) or g W ( Ax ) can be directly extracted. In addition, 
we show how to compute the delayed temporal correla¬ 
tion function, g( 2 \r), from the emission statistics faith¬ 
fully simulating the Hanbury Brown and Twiss setup. 
We explicitly apply the method to non-equilibrium po- 
lariton condensation. 

The model — We consider a multimode bosonic system 
with an energy distribution dictated by the dispersion 
relation E( k) where k = ( k x ,k y ) is the two-dimensional 
wave vector. The system is in contact with a thermal 
reservoir at temperature T and driven by an incoherent 
source with average power P. If only the dominant inter¬ 
actions are considered, the system can be described by 
the Hamiltonian 


T-i = ^kin + ^p-p + Wp-ph + ^pump- (1) 

The first two terms in Eq. 0 describe coherent processes: 
Tt^kin = is the kinetic energy term, in our 

case E k = [Eph (k) — yjE^ h (k) + 4U 2 ]/2 describes the 
lower branch of exciton-polariton dispersion. The pho¬ 
tonic dispersion is given by (k) = ^ 2 /c 2 /(2m p h), the 
exciton-photon Rabi splitting is 2V, and we assume in¬ 
finite exciton mass since m ex m p h- The creation and 
annihilation operators of the bosonic mode with momen¬ 
tum k are denoted by and dk, respectively. The second 
term 


U 


p-p 


b r k 1 k 2 pd kl d ] [ C2 dk 1 +pdk 2 -p (2) 

kik 2 p 


describes the polariton-polariton elastic scattering con¬ 
serving energy and momentum that corresponds to long- 
range interaction. The scattering strength, U klk2P , is 
determined by the excitonic fractions, i.e., the Hop- 
field coefficients, X k , of the initial and final states [45], 
I46J i Lk x k 2 p — b^X kl Xk 2 Xk 1+p Xk 2 — p , where Uq — 
GE^a^/S, Eb and ap are the exciton binding energy and 
the Bohr radius, and S is the system area. 

The last two terms in Eq. 0 describe incoherent pro¬ 
cesses which we treat stochastically. Here, TL P . P h ac¬ 
counts for the interaction of the polaritons with a ther¬ 
mal bath of acoustic phonons. To this end, we introduce 
a Frohlich-type Hamiltonian PHEJ 


u 


p-ph 


E f G 'q«k, «k 2 frq + h.C. 


ki,k 2 


27r 


(3) 


where the sum and integration are performed under 
the condition of energy and momentum conservation, 
|E kl - £ k2 | = do; q , |q| 2 = |ki - k 2 | 2 + q 2 z . The phonons 
are described by the operators b t and 5 q , and their dis¬ 
persion relation, /kj q = huyjq 2 + q 2 + q 2 , is determined 

by the speed of sound u of the material [48]. Above, 
G q is the exciton-phonon interaction strength, the mi¬ 
croscopic derivation of which and typical values can be 
found in Refs. m • We also assume that the phonon 
reservoir remains thermalized, i.e. (5^6 q ) = n p h(7kj q ) = 
{exp[hu/(k B T)} - l} -1 . 

The effect of incoherent pumping is described in the 
rotating-wave approximation by the last term 

^pump = ^ y ^ ^k^dkd^ + f7k£^k^£^ ’ (^) 

k£ 

of Eq.Jll, where d £ and d J are the operators correspond¬ 
ing to the bosonic pumping reservoir in question at tem¬ 
perature Tp and (d^d^) = rip (i£ k ) = {exp[£?k/(^s^p)] — 
1} _1 . The parameters g^ describe the typical linear 
coupling strengths between the system modes and the 
reservoir modes. Further, we assume that g^ = 7k- 
The Hamiltonian Q allows for particle loss which occurs 
with the rate 7k- The losses coming from 0 are mainly 
caused by leakage of photons from the cavity since their 
lifetime is much smaller than that of excitons. Thus we 
can put 7k = 1 /t£ hot , where r£ hot is the lifetime of pho¬ 
tons. Without loss of generality, we take into account 
the natural photon decay from the microcavity using the 
Hamiltonian (4) (see [48 for details). 

Using Eqs. 0 and we can derive a Lindblad-type 
master equation and express the full set of associated 
quantum jump operators mi as 


7 + 

= v / 7k^p(^ k )dj c , 

(5) 

A" 

= \/7k [hp{Ek) + l]a k , 

(6) 

77 2 

= \/7 klk2 ^ph(^ki — -Uk 2 )d ki <2 k2 , 

(7) 


= ^7 klk2 [^phC^ki — ^k 2 ) + l]«ki4 2 , 

(8) 


where X kl > X k2 and we denote the phonon-mediated 
scattering rate as 7 k ^ k2 [48]. Equations ^ and (J6| de¬ 
scribe the polariton pumping and decay. The average 
power fed into the polariton system due to interaction 
with the pumping reservoir is described by P = np 7 k . 
Equations 0 and |8| describe transitions between the 
polariton modes mediated by the phonon reservoir. It 
should be noted that processes Q of phonon emission 
remain even at T = 0 K. 

The method .— The quantum dynamics of the system 
is simulated using the Monte Carlo wave function tech¬ 
nique m- The procedure is based on the evolution of the 
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FIG. 1: (color online), (a) Polariton lower dispersion branch 
(grey line) showing the discrete bosonic modes (dots) used in 
the computations, (b) Occupation of the lowest mode, iVo as 
a function of the pump power P, for several temperatures in 
the range (0 — 20) K (see legend), (c), (d): Occupations of the 
modes as functions of their energies at (c) T = 5 K and (d) 
T — 20 K for several pump powers in range (0.5 — 25) X P t h 
(see legend). The dashed red lines are the fits by a thermal 
distribution. 


system wave function through the Schrodinger equation, 

*^17) = 7 i>), (9) 

with the effective non-Hermitian Hamiltonian 

n = H - T E AA + W E (10) 


— W -7+t -T+ - — -f- 1 -7- 

' 9 ^kika'-'kjka 9 t/ k 1 k 2 *^kik 2 - 


kik 2 


kik 2 


The non-Hermitian part in Eq. (10) results in an ap¬ 


parent decay of the norm . We generate a 

random number, 77, initially and evolve the system by 
Eq. ([ 9 ]). The condition < 77 determines if a 

jump operator occurs or not, see [31:. After each jump, 
we normalize the state \^{t)) again and generate a new 
number 77. A single realization, j, of this protocol yields 
a quantum trajectory, \'4>(t))j, with j = 1, 2,..., TV. Em¬ 
ploying an ensemble of trajectories, we can approximate 
the system density matrix as 


p ^ = — - N -ivtoo^’ 


( 11 ) 


where p(t) is the actual density matrix of the system. 
The expectation value of any system observable O can 
be found from 


(d(t)) = Tr dm 


lim 

TV—>-oo 


{ty [ dm ]} ■ 


( 12 ) 


This method not only allows to significantly reduce the 
memory consumption by evolving a ket vector instead of 
a density matrix but it is also ideal for parallelization 
due to independence of the quantum trajectories. In our 
computations, we truncate the Hilbert space to a chosen 
global number of excitations [43] in addition to the usual 
truncation per mode, which allows to drastically reduce 
the dimension of the Hilbert space with negligible loss of 
accuracy. Due to possible qBEC, the maximum number 
of excitations in the lowest-energy mode is taken several 
times larger than for the other states. 

Results and discussion. — The parameters we consider 
correspond to a GaAs-based microcavity having a cylin¬ 
drical symmetry, with the Rabi splitting 2V = 10 meV, 
771 p h = 5 x 10 -5 77io, polariton lifetime of r cz 1 / 7 *. = 20 
ps, E]j = 10 meV, ap = 10 nm, and S = 100 / 1 m 2 . 
The system symmetry and the correlations we compute 
here allow to consider the radial coordinate k r only which 
would not hold anymore if excited states correlation at 
different angles were under the scope. 

Initially the system is prepared in its vacuum state, 
|t/>(0)). Each trajectory is composed of a 500-ps evolution 
which is sufficiently longer than the time scales of the 
processes involved to reach the steady state. For each 
set of parameters, we average the results over N = 5000 
trajectories. 

Figure [T| shows the occupations of the modes, = 
(al.dk) for different pump powers P and temperatures. 
Our sampling of the dispersion relation is shown in 
Fig. [lja) where we fix the energy step A E = 0.33 meV 
between each mode. We approximate the Hopfield coeffi¬ 
cients as Xk = 1 /V% and hence the polariton-polariton 
scattering strength is fixed to Pkik 2 p = Uq/A. Its value 
is adjusted to result in a typical chemical potential of 
UqNq/4 = 1 meV of the lowest-mode if the latter is com¬ 
pletely filled and this, to compensate for the low par¬ 
ticle number (see outlooks section below). The maxi¬ 
mum number of excitations is fixed to = 15 for 

the lowest-energy mode and to 7V max = 5 for the other 
modes. The polariton-phonon scattering strength is set 
to / 17 k = ^70 = 0.05 meV. Finally, we assume that the 
incoherent pump operator in Eq. © is acting only on the 
highest-energy excited states of the system. Around the 
threshold power, P = P t h, Nq exceeds the population 
of the other modes. As shown in Fig. 1 (b), its value 
monotonically increases with P, faster for lower tem¬ 
peratures, due to the polariton-polariton scattering and 
phonon-assisted energy relaxation. The highest-energy 
states in the dispersion are fed by the incoherent pump 
and play the role of bottleneck modes M- Thus the 













4 



FIG. 2: (color online), (a) Lowest-mode second-order tem¬ 
poral coherence function at zero delay, g^ 2 \ 0), as a function 
of the relative pump power, P/P t h, for different phonon tem¬ 
peratures. (c) Finite-delay coherence g^ 2 \r) calculated for 3 
different pump powers at T = 0 K (see legends). 


latter demonstrate a large occupation even for P > P t h 
(not shown) although as one can see, the population of 
the other modes decreases with increasing energy. The 
higher the pump power, the greater fraction of particles 
is observed to reside in the lowest-energy mode. At the 
lowest investigated pump power, P = 0.5 x P t h, the dis¬ 
tribution reads JV& = Nq exp[—^/(fc^T)]. Using this we 
can extract effective polariton temperatures of T = 7.5 
K and T = 23 K for Figs. 1(c) and 1(d), respectively. 

Figure [2] shows our results for the lowest-mode second- 
order temporal coherence, 

(4 ( 0 ) 4 ( T ) «0 (r) a 0 ( 0 )) 

9 y, \T) = - 7 - 2 -> ( 13 ) 

(4 (0)fio (0)) 

where the averages are defined using Eq. ( [T2| . From T = 
5 K, we observe a clear crossover from thermal statistics 
with g^ (0) = 2 to a coherent state for which g^ (0) = 1 
with increasing P, as seen in Fig. [2|a). With decreasing 
temperature, the coherence appears at lower pump power 
as expected. 

To compute the delayed g^ (r) shown in Fig. (2^b), 
we work on the polariton decay statistics recording the 
full Eq. <© -related event history over long 100 ns tra¬ 
jectories. It allows us to build the probability G^ 2 \t) 
of having two polariton decays within a delay r. When 
normalized to the corresponding Poissonian distribution, 
imposed by the mean steady state occupation, we are 
able to reconstruct the correlation function and confirm 
the onset of temporal coherence revealed by q^ (r) ~ 1 
forP>P th @9]. 

Figure [3] shows the ID approximation of the first-order 
spatial coherence function 




(xi ,Xj)= lim 
t—yoo 


4 t {Xi,t))$ (; Xj,t)Y 


(14) 


between two points at positions Xi and Xj in the steady 
state (t —cxo). Here, ^(x,t) = ^jk^ lkx cik(i). We clearly 


observe the onset of long-range spatial coherence at large 
pump powers whereas the coherence decays on short dis¬ 
tances for low powers. Comparison of cases T = OK and 
T = 20 K, reveals the expected loss of spatial coherence 
with increasing temperature. 
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FIG. 3: (color online). Steady-state first-order spatial coher¬ 
ence function g^\xi,Xj) at T = 0 K (a) below threshold, 
P — 0.5 X P t h, and (b) above threshold, P = 20 X P t h, as 
a function of distance along the sample. Note that periodic 
boundary conditions are imposed by the Fourier transform. 
Panels (c) and (d) show g^(Ax) — g^\ Ax, 0) for various 
pump powers at (c) T = 0 K and (d) T = 20 K, where Ax is 
the distance from the centre of the system. 

Outlook. — In summary, using a stochastic wave func¬ 
tion approach, we have analyzed the quantum properties 
of a non-equilibrium condensate as a function of pump 
intensity and temperature. Our results exhibit all the 
characteristic features associated with the Bose-Einstein 
condensation of such incoherently driven bosonic parti¬ 
cles in contact with a phonon bath. To account in future 
for larger number of bosons and modes, our results can be 
extended by separating the classical field of each mode, 
evolving according to Langevin equations, from the quan¬ 
tum fluctuations that would be treated through the quan¬ 
tum jumps approach with the requirement of a very small 
number of quanta per mode. This will be addressed in a 
separate study. Finally the impact of decoherence in the 
form of pure dephasing can be straightforwardly added 
as a new set of quantum jumps operators m- 
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